R Notebook

setwd("/Users/efraimshine/Desktop/rna_seq/TCGA/")

Warning: The working directory was changed to /Users/efraimshine/Desktop/rna_seq/TCGA inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

data("XenaData")
write.csv(XenaData, "00_tblXenaHubInfo.csv")

## step 5-a
GeneExpectedCnt_toil = XenaGenerate(subset = XenaHostNames == "toilHub") %>%
  XenaFilter(filterCohorts = "TCGA TARGET GTEx") %>%
  XenaFilter(filterDatasets = "TcgaTargetGtex_gene_expected_count");
XenaQuery(GeneExpectedCnt_toil) %>%
  XenaDownload(destdir = "./")

This will check url status, please be patient.
All downloaded files will under directory ./.
The 'trans_slash' option is FALSE, keep same directory structure as Xena.
Creating directories for datasets...
.//TcgaTargetGtex_gene_expected_count.gz, the file has been download!

library("UCSCXenaTools")

=========================================================================================
UCSCXenaTools version 1.4.8
Project URL: https://github.com/ropensci/UCSCXenaTools
Usages: https://cran.r-project.org/web/packages/UCSCXenaTools/vignettes/USCSXenaTools.html

If you use it in published research, please cite:
Wang et al., (2019). The UCSCXenaTools R package: a toolkit for accessing genomics data
  from UCSC Xena platform, from cancer multi-omics to single-cell RNA-seq.
  Journal of Open Source Software, 4(40), 1627, https://doi.org/10.21105/joss.01627
=========================================================================================
                              --Enjoy it--

library("R.utils")

Loading required package: R.oo
Loading required package: R.methodsS3
R.methodsS3 v1.8.2 (2022-06-13 22:00:14 UTC) successfully loaded. See ?R.methodsS3 for help.
R.oo v1.27.0 (2024-11-01 18:00:02 UTC) successfully loaded. See ?R.oo for help.

Attaching package: ‘R.oo’

The following object is masked from ‘package:R.methodsS3’:

    throw

The following objects are masked from ‘package:methods’:

    getClasses, getMethods

The following objects are masked from ‘package:base’:

    attach, detach, load, save

R.utils v2.12.3 (2023-11-18 01:00:02 UTC) successfully loaded. See ?R.utils for help.

Attaching package: ‘R.utils’

The following object is masked from ‘package:utils’:

    timestamp

The following objects are masked from ‘package:base’:

    cat, commandArgs, getOption, isOpen, nullfile, parse, use, warnings

library("dplyr")

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library("data.table")

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
data.table 1.16.4 using 4 threads (see ?getDTthreads).  Latest news: r-datatable.com

Attaching package: ‘data.table’

The following objects are masked from ‘package:dplyr’:

    between, first, last

library("readr")
library("edgeR")

Loading required package: limma

setwd("/Users/efraimshine/Desktop/rna_seq/TCGA")

Warning: The working directory was changed to /Users/efraimshine/Desktop/rna_seq/TCGA inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

samples_gene_expression = fread("TcgaTargetGtex_gene_expected_count.gz");

|--------------------------------------------------|
|==================================================|
Warning: Discarded single-line footer: <<ENSG00000171757.15 5.8056  8.2322  7.2623  4.5850  3.9165  7.3622  7.2790  7.6363  6.7690  6.7673  4.4019  6.5567  5.1735  5.9430  5.1526  6.8088  6.3393  6.0659  8.9498  5.6185  5.4757  6.3232  8.1606  8.5839  5.5689  7.5140  9.1020  5.2593  6.4981  5.4886  6.2024  6.1106  6.6073  6.4373  3.1699  6.0870  6.9506  4.4182  9.3345  7.4155  8.0521  4.5503  6.5664  4.7214  5.5928  6.5469  5.5680  6.3293  4.3125  5.3377  8.5245  4.6809  4.4770  7.9141  6.9937  5.4621  6.8374  4.4790  6.7127  7.8719  7.0444  3.0000  4.9411  4.8354  5.4835  8.5906  8.9245  1.6508  7.761>>

setwd("/Users/efraimshine/Desktop/rna_seq/TCGA")

Warning: The working directory was changed to /Users/efraimshine/Desktop/rna_seq/TCGA inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

create_back_transformed_df <- function(sample_file, gene_expression_data) {
  # Read the sample list
  list_data <- read.csv(sample_file, header = FALSE)
  vector_data <- list_data$V1

  # Filter out excluded samples
  filtered_vector <- vector_data[vector_data %in% colnames(gene_expression_data)]


  # Filter the gene expression data
  filtered_df <- gene_expression_data %>% select("sample", all_of(filtered_vector))
  rownames(filtered_df) <- filtered_df$sample
  
  # Backtransform the numeric columns
  backtransformed_df <- filtered_df %>%
    mutate(across(where(is.numeric), ~ 2^. - 1))
  rownames(backtransformed_df) <- backtransformed_df$sample
  # Set row names and remove the first column

  return(backtransformed_df)
}

g6pd_deleted_samples_file <- "g6pd_deleted_samples.csv"
#sample_file_g6pd_control_normal_expression_file = "g6pd_normal_samples_control_group.csv"
sample_file_g6pd_control_normal_expression_file = "G6PD_normal_cancer_type_based_selection.csv"

# Call the function and store the result
g6pd_del_df <- create_back_transformed_df(g6pd_deleted_samples_file, samples_gene_expression)
g6pd_normal_control_df <- create_back_transformed_df(sample_file_g6pd_control_normal_expression_file, samples_gene_expression)
head(g6pd_normal_control_df[, 1:10], 10)

head(g6pd_del_df[, 1:10], 10)

g6pd_normal_control_df

g6pd_del_df

NA

# checking control data frame samples
# ctr_rnames <- colnames(g6pd_normal_control_df[,-1])
# ctr_orig_df <- read.csv("/Users/efraimshine/Desktop/rna_seq/TCGA/g6pd_high_exp_1.75to2.25_z_score_samples_names.csv")
# ctr_orig_df
# ctr_rnames
# common_names <- intersect(ctr_rnames,ctr_orig_df$Sample.ID)
# common_names
# print(common_names)
g6pd_normal_control_df

# Define the arguments
library("data.table")
library("readr")
library("dplyr")
library("edgeR")




# Add prefixes to column names
colnames(g6pd_del_df) <- paste0("DELE_", colnames(g6pd_del_df))
colnames(g6pd_normal_control_df) <- paste0("NORM_", colnames(g6pd_normal_control_df))

# Combine the data frames by rows
merged_df <- merge(g6pd_del_df, g6pd_normal_control_df, by.x="DELE_sample",by.y = "NORM_sample")
colnames(merged_df)[colnames(merged_df) == "DELE_sample"] <- "sample"
# rownames(merged_df) <- merged_df$sample
# merged_df$sample <- NULL

# rownames(g6pd_del_df) <- g6pd_del_df$DELE_sample
# g6pd_del_df$DELE_sample <- NULL

# rownames(g6pd_normal_control_df) <- g6pd_normal_control_df$NORM_sample
# g6pd_normal_control_df$NORM_sample <- NULL
merged_df

NA
NA

merged_df_rownames <- merged_df[,-1]
rownames(merged_df_rownames) <- merged_df$sample
# merged_df_round_vector <- lapply(merged_df_rownames, round , digits = 2)
# merged_df_rownames_rounded <- data.frame(merged_df_round_vector)
#rownames(merged_df_rownames_rounded) <- rownames(merged_df_rownames)
merged_df_rownames_rounded <- merged_df_rownames %>% mutate(across(where(is.numeric), round, digits = 0))

Warning: There was 1 warning in `mutate()`.
ℹ In argument: `across(where(is.numeric), round, digits = 0)`.
Caused by warning:
! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
Supply arguments directly to `.fns` through an anonymous function instead.

  # Previously
  across(a:b, mean, na.rm = TRUE)

  # Now
  across(a:b, \(x) mean(x, na.rm = TRUE))
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

matrix_for_deseq <- as.matrix(merged_df_rownames_rounded)
rownames(matrix_for_deseq) <- rownames(merged_df_rownames_rounded)
merged_df_rownames

merged_df_rownames_rounded

setwd("/Users/efraimshine/Desktop/rna_seq/TCGA")

Warning: The working directory was changed to /Users/efraimshine/Desktop/rna_seq/TCGA inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

metadata <- data.frame(sample_names <- colnames(merged_df[,-1]))
#row.names(metadata) <- colnames(merged_df)
metadata <- mutate(metadata, G6PD_STATUS = ifelse(row_number() <= length(g6pd_del_df[,-1]), "del", "normal"))

metadata_rnames <- metadata
rownames(metadata_rnames) <- metadata$sample_names....colnames.merged_df
metadata_rnames$sample_names....colnames.merged_df....1.. <- NULL
metadata

metadata_rnames

NA
NA
NA

# metadata <- mutate(metadata, G6PD_STATUS = ifelse(row_number() <= length(g6pd_del_df[,-1]), "del", "normal"))

setwd("/Users/efraimshine/Desktop/rna_seq/TCGA")

Warning: The working directory was changed to /Users/efraimshine/Desktop/rna_seq/TCGA inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

metadata_rnames_for_info_merge <- metadata_rnames
metadata_rnames_for_info_merge

#expanding meta data information for analysis of g6pd deleted samples VS 200 g6pd normal samples from the center of g6pd expression information
library("dplyr")

#samples_info_deleted_vs_200_patients <- read_tsv("g6pd_samples_info_deleted_vs_two_hundreds.tsv") 
#samples_info_deleted_vs_200_patients <- as.data.frame(samples_info_deleted_vs_200_patients) %>% dplyr::select(`Sample ID`,`Cancer Type`,`Cancer Type Detailed`)

samples_info_deleted_vs_cancer_type_based_selection <- read_tsv("info_g6pd_deleted_vs_normal_cancer_type_based_selection.tsv") 

Rows: 1539 Columns: 63── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (43): Study ID, Patient ID, Sample ID, Neoplasm Disease Stage American Joint Committee on Cancer Code, American Joint Comm...
dbl (20): Diagnosis Age, Aneuploidy Score, Last Communication Contact from Initial Pathologic Diagnosis Date, Birth from Initi...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

samples_info_deleted_vs_cancer_type_based_selection <- as.data.frame(samples_info_deleted_vs_cancer_type_based_selection) %>% dplyr::select(`Sample ID`,`Cancer Type`,`Cancer Type Detailed`,`Study ID`)


# merege metadata with info df
row.names(metadata_rnames_for_info_merge) <- gsub("^[^_]*_", "", row.names(metadata_rnames_for_info_merge))

metadata_rnames_merged_with_info <- merge(metadata_rnames_for_info_merge,samples_info_deleted_vs_cancer_type_based_selection, by.x = "row.names", by.y = "Sample ID",all.x = TRUE)
row.names(metadata_rnames_merged_with_info) <- metadata_rnames_merged_with_info$Row.names
metadata_rnames_merged_with_info$Row.names <- NULL
metadata_rnames_merged_with_info <- metadata_rnames_merged_with_info[order(metadata_rnames_merged_with_info$G6PD_STATUS), ]
metadata_rnames_merged_with_info_orderd <- metadata_rnames_merged_with_info[match(row.names(metadata_rnames_for_info_merge), row.names(metadata_rnames_merged_with_info)), ]

row.names(metadata_rnames_merged_with_info_orderd) <- row.names(metadata_rnames)
colnames(metadata_rnames_merged_with_info_orderd)[2] <- "cancer_type"
colnames(metadata_rnames_merged_with_info_orderd)[3] <- "cancer_type_detailed"

metadata_rnames_merged_with_info_orderd$cancer_type <- as.factor(metadata_rnames_merged_with_info_orderd$cancer_type)
metadata_rnames_merged_with_info_orderd$cancer_type_detailed <- as.factor(metadata_rnames_merged_with_info_orderd$cancer_type_detailed)
metadata_rnames_merged_with_info

metadata_rnames_merged_with_info_orderd

NA

colnames(metadata_rnames_merged_with_info_orderd)[4] <- "study_ID"
metadata_rnames_merged_with_info_orderd$cancer_type <- as.factor(metadata_rnames_merged_with_info_orderd$cancer_type)
metadata_rnames_merged_with_info_orderd$cancer_type_detailed <- as.factor(metadata_rnames_merged_with_info_orderd$cancer_type_detailed)
metadata_rnames_merged_with_info_orderd$G6PD_STATUS <- as.factor(metadata_rnames_merged_with_info_orderd$G6PD_STATUS)
metadata_rnames_merged_with_info_orderd$`Study ID` <- as.factor(metadata_rnames_merged_with_info_orderd$`study_ID`)
metadata_rnames_merged_with_info

metadata_rnames_merged_with_info_orderd

names_df <- as.data.frame(data_col <- colnames(matrix_for_deseq),metadatarows <- rownames(metadata_rnames))
names_df$`data_col <- colnames(matrix_for_deseq)` == rownames(names_df)

   [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
  [25] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
  [49] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
  [73] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
  [97] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [121] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [145] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [169] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [193] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [217] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [241] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [265] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [289] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [313] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [337] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [361] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [385] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [409] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [433] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [457] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [481] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [505] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [529] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [553] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [577] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [601] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [625] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [649] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [673] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [697] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [721] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [745] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [769] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [793] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [817] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [841] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [865] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [889] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [913] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [937] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [961] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [985] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [ reached getOption("max.print") -- omitted 428 entries ]

names_df

unique_cancers <- unique(metadata_rnames_merged_with_info_orderd$cancer_type)

results_list <- list()

for (cancer in unique_cancers) {
  subset_meta <- metadata_rnames_merged_with_info_orderd[metadata_rnames_merged_with_info_orderd$cancer_type == cancer, ]
  subset_counts <- matrix_for_deseq[, rownames(subset_meta)]
  
  dds_loop <- DESeqDataSetFromMatrix(countData = subset_counts,
                                colData = subset_meta,
                                design = ~ G6PD_STATUS)
  dds_loop <- DESeq(dds_loop)
  
  results_list[[cancer]] <- results(dds_loop, contrast = c("G6PD_STATUS", "del", "normal"))
}

converting counts to integer mode
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 18811 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
converting counts to integer mode
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 22176 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
converting counts to integer mode
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 3719 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
converting counts to integer mode
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 22515 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
converting counts to integer mode
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
converting counts to integer mode
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 3271 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
converting counts to integer mode
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 19251 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
converting counts to integer mode
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 2076 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
converting counts to integer mode
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 2299 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
converting counts to integer mode
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 2671 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
converting counts to integer mode
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 2497 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
converting counts to integer mode
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 23669 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
converting counts to integer mode
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 3018 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
converting counts to integer mode
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 1024 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
converting counts to integer mode
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 23071 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
converting counts to integer mode
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 23429 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing

# Load required libraries
library(DESeq2)
library(clusterProfiler)
library(org.Hs.eg.db)
library(ggplot2)

# Initialize lists to store results
filtered_results_list <- list()
GO_enrichment_results <- list()

# Loop through each cancer type in results_list
for (cancer in names(results_list)) {
  
  # Extract results (already a DESeqResults object)
  res05 <- results_list[[cancer]]
  
  # Summary of DEGs
  cat("\nSummary for", cancer, ":\n")
  print(summary(res05))
  
  # Generate MA plot
  png(paste0("MA_plot_", cancer, ".png"), width = 800, height = 600)
  plotMA(res05, ylim=c(-2,2), main=paste("MA Plot -", cancer))
  dev.off()
  
  # Remove NA values
  sigs <- na.omit(res05)
  
  # Filter significantly differentially expressed genes (padj < 0.05 & baseMean > 50)
  sigs <- sigs[sigs$padj < 0.05 & sigs$baseMean > 0,]
  
  # Store filtered results
  filtered_results_list[[cancer]] <- sigs
  
  # Extract genes with log2FoldChange > 0.5
  genes_to_test <- rownames(sigs[sigs$log2FoldChange > 0,])
  
  # Remove Ensembl version numbers (if present)
  genes_to_test_no_dots <- gsub("\\..*", "", genes_to_test)
  
  # Perform GO enrichment analysis (if there are genes to test)
  if (length(genes_to_test_no_dots) > 0) {
    GO_results <- enrichGO(gene = genes_to_test_no_dots, 
                           OrgDb = "org.Hs.eg.db", 
                           keyType = "ENSEMBL", 
                           ont = "BP", 
                           pAdjustMethod = "BH",
                           pvalueCutoff = 0.05)
    
    # Convert to dataframe and store results
    GO_enrichment_results[[cancer]] <- as.data.frame(GO_results)
  }
}

Summary for Adrenocortical Carcinoma :

out of 24498 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 23, 0.094%
LFC < 0 (down)     : 118, 0.48%
outliers [1]       : 0, 0%
low counts [2]     : 12484, 51%
(mean count < 12)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

NULL

Summary for Bladder Cancer :

out of 29831 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 55, 0.18%
LFC < 0 (down)     : 338, 1.1%
outliers [1]       : 0, 0%
low counts [2]     : 17672, 59%
(mean count < 18)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

NULL

Summary for Breast Cancer :

out of 32003 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 23, 0.072%
LFC < 0 (down)     : 32, 0.1%
outliers [1]       : 381, 1.2%
low counts [2]     : 13474, 42%
(mean count < 1)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

NULL

Summary for Colorectal Cancer :

out of 28413 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 199, 0.7%
LFC < 0 (down)     : 294, 1%
outliers [1]       : 0, 0%
low counts [2]     : 16250, 57%
(mean count < 13)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

NULL

Summary for Mature B-Cell Neoplasms :

out of 22424 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 156, 0.7%
LFC < 0 (down)     : 151, 0.67%
outliers [1]       : 488, 2.2%
low counts [2]     : 8767, 39%
(mean count < 4)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

NULL

Summary for Esophagogastric Cancer :

out of 33539 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 18, 0.054%
LFC < 0 (down)     : 24, 0.072%
outliers [1]       : 664, 2%
low counts [2]     : 17135, 51%
(mean count < 5)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

NULL

Summary for Glioblastoma :

out of 27196 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 191, 0.7%
LFC < 0 (down)     : 274, 1%
outliers [1]       : 0, 0%
low counts [2]     : 12906, 47%
(mean count < 5)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

NULL

Summary for Head and Neck Cancer :

out of 30595 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 1136, 3.7%
LFC < 0 (down)     : 2044, 6.7%
outliers [1]       : 0, 0%
low counts [2]     : 15811, 52%
(mean count < 3)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

NULL

Summary for Hepatobiliary Cancer :

out of 27633 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 2, 0.0072%
LFC < 0 (down)     : 40, 0.14%
outliers [1]       : 331, 1.2%
low counts [2]     : 11533, 42%
(mean count < 1)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

NULL

Summary for Non-Small Cell Lung Cancer :

out of 30237 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 24, 0.079%
LFC < 0 (down)     : 25, 0.083%
outliers [1]       : 614, 2%
low counts [2]     : 13208, 44%
(mean count < 2)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

NULL

Summary for Ovarian Epithelial Tumor :

out of 30526 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 269, 0.88%
LFC < 0 (down)     : 332, 1.1%
outliers [1]       : 0, 0%
low counts [2]     : 14019, 46%
(mean count < 4)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

NULL

Summary for Prostate Cancer :

out of 29672 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 8, 0.027%
LFC < 0 (down)     : 30, 0.1%
outliers [1]       : 0, 0%
low counts [2]     : 14182, 48%
(mean count < 3)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

NULL

Summary for Sarcoma :

out of 28910 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 162, 0.56%
LFC < 0 (down)     : 141, 0.49%
outliers [1]       : 1152, 4%
low counts [2]     : 11467, 40%
(mean count < 2)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

NULL

Summary for Thyroid Cancer :

out of 29102 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 11, 0.038%
LFC < 0 (down)     : 6, 0.021%
outliers [1]       : 180, 0.62%
low counts [2]     : 11631, 40%
(mean count < 1)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

NULL

Summary for Endometrial Cancer :

out of 29480 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 22, 0.075%
LFC < 0 (down)     : 53, 0.18%
outliers [1]       : 0, 0%
low counts [2]     : 15315, 52%
(mean count < 7)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

NULL

Summary for Melanoma :

out of 29825 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 4, 0.013%
LFC < 0 (down)     : 83, 0.28%
outliers [1]       : 0, 0%
low counts [2]     : 15945, 53%
(mean count < 7)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

NULL

# Save results
save(filtered_results_list, GO_enrichment_results, file = "DEG_GO_results.RData")

# Print message when done
cat("Analysis complete! Filtered DEGs and GO enrichment results are saved.\n")

Analysis complete! Filtered DEGs and GO enrichment results are saved.

multi_results <- load("DEG_GO_results.RData")
GO_enrichment_results[["Mature B-Cell Neoplasms"]]

library("DESeq2")
dds_multi <- DESeqDataSetFromMatrix(countData = matrix_for_deseq,
                              colData = metadata_rnames_merged_with_info_orderd,
                              design = ~ cancer_type + G6PD_STATUS + cancer_type:G6PD_STATUS)

converting counts to integer mode
  Note: levels of factors in the design contain characters other than
  letters, numbers, '_' and '.'. It is recommended (but not required) to use
  only letters, numbers, and delimiters '_' or '.', as these are safe characters
  for column names in R. [This is a message, not a warning or an error]

dds_multi$G6PD_STATUS <- relevel(dds_multi$G6PD_STATUS, ref = "normal")
dds_multi

class: DESeqDataSet 
dim: 38608 1428 
metadata(1): version
assays(1): counts
rownames(38608): ENSG00000000938.12 ENSG00000000971.15 ... ENSGR0000280767.2 ENSGR0000281849.2
rowData names(0):
colnames(1428): DELE_TCGA-OR-A5LC-01 DELE_TCGA-DK-A1AE-01 ... NORM_TCGA-ZS-A9CD-01 NORM_TCGA-ZS-A9CF-01
colData names(5): G6PD_STATUS cancer_type cancer_type_detailed study_ID Study ID

dds_multi <- DESeq(dds_multi)

estimating size factors
  Note: levels of factors in the design contain characters other than
  letters, numbers, '_' and '.'. It is recommended (but not required) to use
  only letters, numbers, and delimiters '_' or '.', as these are safe characters
  for column names in R. [This is a message, not a warning or an error]
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
  Note: levels of factors in the design contain characters other than
  letters, numbers, '_' and '.'. It is recommended (but not required) to use
  only letters, numbers, and delimiters '_' or '.', as these are safe characters
  for column names in R. [This is a message, not a warning or an error]
final dispersion estimates
fitting model and testing
1615 rows did not converge in beta, labelled in mcols(object)$betaConv. Use larger maxit argument with nbinomWaldTest
  Note: levels of factors in the design contain characters other than
  letters, numbers, '_' and '.'. It is recommended (but not required) to use
  only letters, numbers, and delimiters '_' or '.', as these are safe characters
  for column names in R. [This is a message, not a warning or an error]
-- replacing outliers and refitting for 22497 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
Error: vector memory limit of 16.0 Gb reached, see mem.maxVSize()

resultsNames(dds)

resOrdered <- res[order(res$pvalue),]
summary(res)

sum(res$padj < 0.05, na.rm=TRUE)

res05 <- results(dds, alpha=0.05)
summary(res05)

plotMA(res, ylim=c(-2,2))

resSig <- subset(resOrdered, padj < 0.05)
resSig

#filtering the results
results_df <- read.csv("/Users/efraimshine/Desktop/rna_seq/TCGA/G6PD_abovepadj_0.05.csv")
results_df
results_filtered <- results_df %>% filter(log2FoldChange >0 ,padj < 0.05)
results_filtered
write.csv(as.data.frame(results_filtered), 
          file="/Users/efraimshine/Desktop/rna_seq/TCGA/deseq2_res_g6pd/G6PD_results_filtered.csv")

# clustering
library(clusterProfiler)
library(org.Hs.eg.db)
library(AnnotationDbi)

# results as in tutorial notebook
sigs <- na.omit(res)
sigs
sigs <- sigs[sigs$padj < 0.05 & sigs$baseMean > 50,]
sigs
genes_to_test <- rownames(sigs[sigs$log2FoldChange > 0.5,]) 

genes_to_test_no_dots <- gsub("\\..*", "", genes_to_test)

#GO analysis
GO_results <- enrichGO(gene = genes_to_test_no_dots, OrgDb = "org.Hs.eg.db", keyType = "ENSEMBL", ont = "BP")
as.data.frame(GO_results)

fit <- plot(barplot(GO_results, showCategory = 20))
png("out.png", res = 250, width = 1400, height = 1800)
print(fit)
dev.off()

# GO_results <- enricher()
# as.data.frame(GO_results)

# # prepaer meta data for multi factor analysis 
# metadata_rnames_merged_with_info_orderd
# metadata_rnames_merged_with_info_orderd$G6PD_STATUS <- factor(metadata_rnames_merged_with_info_orderd$G6PD_STATUS)
# table(metadata_rnames_merged_with_info_orderd$cancer_type, metadata_rnames_merged_with_info_orderd$G6PD_STATUS)
# metadata_rnames_merged_with_info_orderd$cancer_type <- droplevels(metadata_rnames_merged_with_info_orderd$cancer_type)
# metadata_rnames_merged_with_info_orderd$G6PD_STATUS <- droplevels(metadata_rnames_merged_with_info_orderd$G6PD_STATUS)
# metadata_rnames_merged_with_info_orderd

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.